DNA Metabarcoding Highlights Cyanobacteria as the Main Source of Primary Production in a Pelagic Food Web Model - Supplementary information: Model analysis
Andreas Novotny, Baptiste Serandour, Susanne Kortsch, Benoit Gauzens, Kinlan Mehdi Goulwen Jan, and Monika Winder
Contact: Andreas Novotny, mail@andreasnovotny.se
library(tidyverse)
library(xts)
library(igraph)
library(svglite)
library(circlize)
library(lubridate)
library(forcats)
library(fluxweb)
library(car)1 Data Preparations
Figure: Overview of data preparation methods
1.1 Selectivity index based on DNA metabarcoding
Individuals of zooplankton were sampled at three locations in the Baltic Sea (Landsort Deep, Gotland Deep, and Bornholm Deep) over the season between 2017 and 2020 (Fig. E1, Tab. E3). Zooplankton samples were collected with vertical hauls from 0-30 and 30-60 m using a 90 µm-WP2 plankton net. For estimating available prey composition, water samples were collected with 10L Niskin bottles with 5 m depth intervals above the thermocline of 30 m depth or with a 20 m long hose. The depths were mixed by adding an equal volume of water from the Niskin bottles. 1-3L were sequentially filtered onto 25 mm diameter filters with 20 µm, 2 µm (polycarbonate), and 0.2 µm (nylon) pore sizes. Sampling for DNA analysis followed the protocol of the Swedish national monitoring program for microscopic count data.
We used a metabarcoding assay of the V3-V4 region of 16S rRNA gene of DNA extracted from the dominant zooplankton species, including Acartia spp., Centropages hamatus, Evadne nordmanni, Pseudocalanus spp., Temora longicornis, Synchaeta baltica, Bosmina spp., Keratella spp., and the water samples. Libraries were constructed using primers 341F-805R (Herlemann et al. 2011;
et al. 2016) and the sequences were assigned to the SILVA (Pruesse et al. 2007) and PhytoREF databases (Decelle et al. 2015).
Load DNA dataset zooplankton-consumer and prey availability data:
zoopl_16S <- readRDS("./Data/16SrRNA_zp.rds") %>%
dplyr::mutate(StationDate = paste(STATION_ID, SAMPLE_date),
StationDateOTU = paste(STATION_ID, SAMPLE_date, OTU),
NODE_NAME = ifelse(SORTED_genus == "S.baltica", "Synchaeta", SORTED_genus))
head(zoopl_16S)Create table of samples, to get an overview of the number of sample replicates per day, station, and zooplankton species.
zoopl_16S %>%
filter(SORTED_genus != "Water") %>%
group_by(Sample, SORTED_genus, SAMPLE_date, STATION_ID) %>%
summarize() %>%
group_by(SORTED_genus, SAMPLE_date, STATION_ID) %>%
summarize(n = length(Sample)) %>%
write_csv("./Outputs/Tab.S3.csv")## `summarise()` has grouped output by 'Sample', 'SORTED_genus', 'SAMPLE_date'.
## You can override using the `.groups` argument.
## `summarise()` has grouped output by 'SORTED_genus', 'SAMPLE_date'. You can
## override using the `.groups` argument.To assess the selectivity index S(i) for each prey taxon (i) of a consumer, we use the standardized forage ratio as described by (Chesson 1983). S(i) = Rg(i)/Re(i)/sum:k(Rg(k)/Re(k)) Rg(i) is the Read count of i (Abundance) in the consumers gut . Re(i) is the Read count of i in the environmental sample.
Create a paired data set, with each OTU in both the water and gut sample on one row.
environment <- zoopl_16S %>%
# Select Water samples only
dplyr::filter(NODE_NAME == "Water")
paired_data <- zoopl_16S %>%
# Select gut samples only
dplyr::filter(NODE_NAME != "Water") %>%
# Merge the Zooplankton and the water dataset, on the OTU level.
dplyr::left_join(environment, by="OTU") %>%
# Remove OTUs with no read in the water sample.
dplyr::filter(is.na(Sample.y)==FALSE,
Abundance.y>0)
head(paired_data)Calculate Selectivity indices for each OTU and sample
selectivity_data <- paired_data %>%
# Calculate preference value based on water Rg(i)/Re(i)
dplyr::mutate(Preference = Abundance.x/Abundance.y) %>%
# Summarize per sample Preference/sum(prefeerence) ; RA(ij)/RA(iw) / sum(k, RA(kj)/RA(kw))
dplyr::group_by(Sample.x) %>%
dplyr::summarize(Selectivity_index=Preference/sum(Preference),
OTU = OTU, NODE_NAME = NODE_NAME.x,
Domain=Domain.x, Supergroup=Supergroup.x, Class = Class.x,
Order=Order.x, Family=Family.x, Genus=Genus.x, Species=Species.x, SAMPLE_date=SAMPLE_date.x)
head(selectivity_data)Summarize Selectivity indices for sample replicates and OTUs.
summarized_selectivity <- selectivity_data %>%
# Mean Si for all OTUs for sample replicates.
dplyr::group_by(NODE_NAME, Class, Order, Family, Genus, Species, OTU, SAMPLE_date) %>%
dplyr::summarise(Selectivity_index= mean(Selectivity_index)) %>%
# Merge OTUs at suitable level. (Here, Class)
dplyr::group_by(NODE_NAME, Order, SAMPLE_date) %>%
dplyr::summarise(Selectivity_index= sum(Selectivity_index)) %>%
#Filter uniportant contributions
dplyr::filter(Selectivity_index > 0.00) %>%
dplyr::group_by(NODE_NAME, SAMPLE_date) %>%
dplyr::summarize(Selectivity_index=Selectivity_index/sum(Selectivity_index), Order)
head(summarized_selectivity)rm(environment, paired_data, zoopl_16S)1.2 Biomasses
The population biomasses of zooplankton and phytoplankton were retrieved from the Swedish national pelagic phytoplankton and zooplankton monitoring with sampling intensity ranging between monthly and weekly samples(51). Individual body masses of phytoplankton and zooplankton were retrieved from COMBINE guidelines for plankton monitoring in the Baltic Sea(52). We calculated daily biomass estimates over one year by linearly interpolating data from samples taken between 2007 and 2018.
Zooplankton <- readRDS("Data/Zooplankton_shark.rds")
Phytoplankton <- readRDS("./Data/Phytoplankton_shark.rds")Daily interpolation of biomass
We calculated daily biomass estimates by linearly interpolating data from samples taken between 2007 and 2018. First we define the function for daily interpolation:
dailyInterpretation <- function(data, taxa="Genus") {
require(tidyverse, xts)
allDates <- seq.Date(
#Creating a sequence of days for dataset
min(as.Date(data$SDATE)),
max(as.Date(data$SDATE)),
"day")
daydata <- data %>%
# Turn into "matrix" format
spread(key={{taxa}}, value=Value) %>%
# Replaceing Nas with 0 for actual samples
replace(is.na(.), 0) %>%
# Adding the sequence of days to existing dataset
full_join(x=data.frame(SDATE=allDates), y=.) %>%
# Order acc to date
arrange(SDATE) %>%
# Linnear na approximation of all numeric variables in dataset
mutate_if(is.numeric, funs(na.approx(.))) %>%
# Back to long format
gather(key = "Taxa", colnames(.)[-1], value = "Value") %>%
# Summarize days of each year
mutate(Abundance = as.numeric(as.character(Value)),
DOY = as.numeric(strftime(SDATE, format = "%j")),
Year = as.numeric(strftime(SDATE, format = "%Y"))) %>%
group_by(DOY, Taxa) %>%
summarize(Abundance = mean(Abundance))
return(daydata)
}Daily interpolation of zooplankton biomass
dip_zooplankton <- function(x) {
zoopl <- Zooplankton %>%
dplyr::filter(Station == x,
Depth %in% c(30,60),
Year > 2006,
Parameter == "Wet weight/volume",
dev_stage_code != "NP",
Genus %in% c(
"Evadne",
"Bosmina",
"Acartia",
"Temora",
"Centropages",
"Keratella",
"Synchaeta",
"Pseudocalanus"
)) %>%
group_by(SDATE, Genus) %>%
summarize(Value=sum(Value)*30/1000) %>% # *30 To get the full water column, /1000 to get from mg to g.
dailyInterpretation(taxa="Genus") %>%
mutate(NODE_NAME=Taxa,
Station=x)
return(zoopl)
}
zooplankton_dip <- dip_zooplankton("BY31 LANDSORTSDJ") %>%
bind_rows(dip_zooplankton("BY15 GOTLANDSDJ")) %>%
bind_rows(dip_zooplankton("BY5 BORNHOLMSDJ"))## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
##
## # Simple named list: list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
##
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))head(zooplankton_dip)Daily interpolation of phytoplankton
dip_phytoplankton <- function(x) {
phytopl <- Phytoplankton %>%
filter(Station == x,
Year > 2006,
unit == "ugC/l") %>%
group_by(SDATE, Order, Depth) %>%
summarize(Value=sum(Value)) %>%
group_by(SDATE, Order) %>%
summarize(Value=mean(Value)*0.12) %>%
# *0.12 to correct for depth *30, carbon content*4, and parameter *0.001 ug/l -> g/m
dailyInterpretation(taxa="Order") %>%
mutate(NODE_NAME=Taxa,
Station=x)
}
phytoplankton_dip <- dip_phytoplankton("BY31 LANDSORTSDJ") %>%
bind_rows(dip_phytoplankton("BY15 GOTLANDSDJ")) %>%
bind_rows(dip_phytoplankton("BY5 BORNHOLMSDJ"))## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
## `.name_repair` is omitted as of tibble 2.0.0.
## ℹ Using compatibility `.name_repair`.
## ℹ The deprecated feature was likely used in the tidyr package.
## Please report the issue at <]8;;https://github.com/tidyverse/tidyr/issueshttps://github.com/tidyverse/tidyr/issues]8;;>.head(phytoplankton_dip)1.3 Combine data sets to Links and Edges
In this section we match the different data sources to create a data frame describing the edges (organism groups) of the food web, and one matrix describing the links (feeding selectivity) of the edges.
Filter and merge biomass data
# List prey orders detected in zooplankton guts
zp_prey <- summarized_selectivity %>%
pull(Order) %>%
unique()
# Filter biomass data by detected prey only
phytoplankton_abundance <- phytoplankton_dip %>%
filter(NODE_NAME %in% zp_prey) %>%
select(DOY, Abundance, NODE_NAME, Station)
# Plot all sources per station
phytoplankton_abundance %>%
ggplot(aes(DOY, Abundance)) +
geom_line(aes(color=NODE_NAME)) +
facet_wrap(~Station)
head(phytoplankton_abundance)Biomasses <- zooplankton_dip %>%
select(DOY, Abundance, NODE_NAME, Station) %>%
union(phytoplankton_abundance)Filter selectivity data, Links
# List potential prey that is detected or monitored in the water microscopically
pp <- phytoplankton_abundance %>%
pull(NODE_NAME) %>%
unique()
# Filter selectivity dataset by prey availability
Links <- summarized_selectivity %>%
filter(Order %in% pp) %>%
mutate(Predator = NODE_NAME, Prey = Order) %>%
group_by(Predator, Prey) %>%
summarise(Preference = mean(Selectivity_index))
# Save Selectivity index table
pivot_wider(Links, names_from = Prey, values_from = Preference, values_fill = 0) %>%
relocate(c("Predator",rev(ord.prey.fw))) %>%
arrange(match(Predator, ord.pred.fw)) %>%
write_csv("./Outputs/Tab.S2.csv")
head(Links)1.4 Construct data with information of nodes
bodymasses <- tibble(NODE_NAME = c("Pseudocalanus",
"Temora",
"Centropages",
"Acartia",
"Evadne",
"Bosmina",
"Keratella",
"Synchaeta"),
BODYMASS = c(1.69E-05,
1.69E-05,
1.69E-05,
9.50E-06,
3.44E-05,
3.44E-05,
1.00E-06,
6.00E-06))
nodenames <- Biomasses %>%
pull(NODE_NAME) %>%
unique() %>% tibble(NODE_NAME = .)
Nodes <- merge(bodymasses, nodenames, by=c("NODE_NAME"), all = T) %>%
mutate(BODYMASS = ifelse(NODE_NAME %in% bodymasses$NODE_NAME, as.numeric(BODYMASS), 1),
ORG_TYPE = ifelse(NODE_NAME %in% bodymasses$NODE_NAME, "Animal", "Plant"))
head(Nodes)1.5 Save all modified datasets
saveRDS(Links, "./Data/Mod/Links.rds")
saveRDS(Nodes, "./Data/Mod/Nodes.rds")
saveRDS(Biomasses, "./Data/Mod/Biomasses.rds")2 Model Definition
2.1 Define FoodWebModel Function
The metabolic rate for each node was calculated based on metabolic scaling theory (4), where the metabolic rate Xi (J/s) is derived from species’ body masses Mi and an allometric scaling constant, i.e., a normalization constant a =17.17 for invertebrates and a = 18.47 for vertebrates. Bei represents the total population biomass per area unit (g/m2) and is multiplied by the metabolic rate (flux per gram biomass) to obtain a metabolic rate estimate at the population level. Metabolic rates are also adjusted to ambient temperature, where E is the activation energy, K is Boltzmann’s constant, and T is the absolute temperature in Kelvin. We used a bioenergetic model described by29 to calculate energy fluxes (J/day/m2) between all nodes in the food web for each day of the year, and each of the three stations. The model builds on a steady-state assumption where the energy consumption of each species equals the energetic losses (Gi = Li). In this model, losses of each population (Li) were defined as the sum of metabolic losses for the population, and the energy flux lost to predation. Gains for each population were calculated as the sum of all fluxes from prey to predator multiplied with the assimilation efficiency, i.e., conversion of consumed biomass into energy and an absolute prey preference constant. The efficiency was based on the functional group of the prey and put to 0.7755. The prey preference was based on the selectivity index calculated from DNA abundances scaled with the relative biomass abundance of the prey species.
#' Estimating fluxes of energy
#' For a given day and time, calculates fluxes in the food web,
#' This is a wrapper function for fluxweb::fluxing, not intended for export.
#' @param nodes Nodes dataframe Should contain: BODYMASS, ORG_TYPE, and NODE_NAME
#' @param links DataFrame with links: Predator, Prey, Preference
#' @param biomasses Dataframe with biomass information: DOY, Abundance, NODE_NAME and Station
#' @param temperatures Dataframe with temperature: DOY and Temperature
#' @param doy Integer: What day of the year?
#' @param station Character string: Name of the station
#' @value data.frame
FoodWebModel <- function(nodes=Nodes,
links=Links,
biomasses=Biomasses,
temperatures=Temperatures,
doy=280,
station="BY31 LANDSORTSDJ") {
# 1. Construct two support functions:
#' Support function: Construct Adjacency matrix from edgelist
#' @param links A dataframe with edges containing col1: Pred, Col2: Prey.
#' @param weight Column name that contains whe edge weight.
#' @value numeric.matrix
edgelistToAdjacency <- function(links=Links, weight="Preference") {
graph <- graph.data.frame(links)
mat <- get.adjacency(graph, attr=weight) %>% as.matrix()
adjacency <- mat %>%
unlist() %>%
as.numeric() %>%
matrix(nrow=nrow(mat)) %>%
t()
return(adjacency)
} # End of edgelistToAdjacency()
#' Support function: Convert matrix to long format dataframe
#' @param mat a matrix
#' @value tibble. Col1: rownames(mat), col2: colnames(mat), col3: mat.
matToDf <- function(mat) {
DF <- list(Row = rownames(mat)[row(mat)] %||% row(mat),
Col = colnames(mat)[col(mat)] %||% col(mat),
value = mat) %>%
map_dfc(as.vector) %>%
as_tibble()
return(DF)
} # End of matToDf()
# 2. Calculate LOSSES and define EFFICIENCIES based on metabolic type and temperature
temp <- temperatures %>%
dplyr::filter(DOY==doy) %>%
dplyr::pull(Temperature)
boltz <- 0.00008617343 # Boltzmann constant
tkonst <- 0.69/(boltz*(273.15+temp)) #Temperature metabolic constant
nodes_tmp <- nodes %>%
dplyr::mutate(LOSSES = exp(-0.29 * log(BODYMASS) +
dplyr::recode(ORG_TYPE ,
"Animal" = 17.17,
"ectotherm vertebrates" = 18.47,
"Plant" = 1) -tkonst),
EFFICIENCIES = dplyr::recode(ORG_TYPE,
"Plant" = 0.77,
"Detritus" = 0.4,
"Animal" = 0.906,
"Bacteria" = 0.906))
# Note: Detritus derived from:Gergs, R. and K. O. Rothhaupt (2008) in Freshwater
# Biology 53: 2494-250.However, in practice all resources in this model are "Plants",
# and all resources are "Animals", non ectothermal. All other numbers are not utilized.
# 3. Add BIOMASS to node data frame
nodes_doy <- biomasses %>%
dplyr::filter(DOY==doy, Station==station) %>%
left_join(nodes_tmp, by="NODE_NAME")
# 4. Rearrange dataframe to overlap with LINKS data.
fluxnames <- links %>%
graph.data.frame() %>%
get.adjacency(attr="Preference") %>%
as.matrix() %>%
colnames()
nodes_doy <- nodes_doy[match(fluxnames, nodes_doy$NODE_NAME),] %>%
as.data.frame()
row.names(nodes_doy) <- nodes_doy$NODE_NAME
# 5. Calculate energy fluxes usin the fluxweb packae
flux <- links %>%
edgelistToAdjacency() %>% # Make square numeric matrix.
fluxing(biomasses = nodes_doy$Abundance,
losses = nodes_doy$LOSSES,
efficiencies = nodes_doy$EFFICIENCIES,
bioms.prefs = TRUE,
ef.level = "prey",
bioms.losses = TRUE)
# 6. Clean up the output, and return it to dataframe formate.
colnames(flux) <- fluxnames
row.names(flux) <- fluxnames
food_web <- flux %>%
matToDf() %>%
mutate(DOY=doy,
Station=station,
Predator = Col,
Prey = Row,
Flux = value) %>%
dplyr::select(Predator, Prey, Flux, DOY, Station)
return(food_web)
} # End of FoodWebModel()2.2 Computation and Calculations of metrics
Read in pre filtered data sets that the model is based on.
Links <- readRDS("./Data/Mod/Links.rds")
Nodes <- readRDS("./Data/Mod/Nodes.rds")
Biomasses <- readRDS("./Data/Mod/Biomasses.rds")
Temperatures <- readRDS("./Data/Temperatures.rds")Calculate fluxes for for all days and stations
We calculated fluxes (J/day/m2) between all nodes in the food webs for each day of the year and each station. Here, the above defined FoodWebModel() function is applied to day 30:340, and appended to the three stations, one after the other.
full_model_df <-
mapply(function(Station) {
bind_rows(mapply(FoodWebModel, doy=30:340,
MoreArgs = list(station = Station,
nodes = Nodes,
links = Links,
biomasses = Biomasses,
temperatures = Temperatures),
SIMPLIFY = FALSE))},
Station=c("BY31 LANDSORTSDJ",
"BY15 GOTLANDSDJ",
"BY5 BORNHOLMSDJ"),
SIMPLIFY = FALSE) %>%
bind_rows() %>%
mutate(Flux=Flux*86.4) #Modify unit from j/m2/s to Kj/m2/day
head(full_model_df)Station <- c("BY31 LANDSORTSDJ", "BY15 GOTLANDSDJ", "BY5 BORNHOLMSDJ")
Doy <- 30:340
g <- expand.grid(Station, Doy)
full_model_df <-
mapply(FoodWebModel,
station = g[, 1],
doy = g[, 2],
MoreArgs = list(links=Links,
biomasses=Biomasses,
temperatures=Temperatures),
SIMPLIFY = FALSE) %>%
bind_rows() %>%
mutate(Flux=Flux*86.4) #Modify unit from j/m2/s to Kj/m2/dayCalculate metrics per day
For each node and day of the year, we calculated the total predation losses, normalized predation pressure, and total consumption. We also calculated annual food web metrics by summarizing all daily flux networks (J/year/m2).
Total
# Consumption per predator
Pred_flux <- full_model_df %>%
mutate(NODE_NAME = Predator) %>%
group_by(NODE_NAME, DOY, Station) %>%
summarize(Consumption=sum(Flux))
head(Pred_flux)Summary_statistics_stations <- full_model_df %>%
mutate(NODE_NAME = Prey) %>%
# Calculate predation loss
group_by(NODE_NAME, DOY, Station) %>%
summarize(Predation_loss=sum(Flux)) %>%
left_join(Biomasses) %>%
left_join(Pred_flux) %>%
mutate(Predation_pressure=Predation_loss/Abundance) %>%
gather(key = "Parameter", value = "Value",
Abundance, Predation_loss, Predation_pressure, Consumption)
Summary_statistics_all <- Summary_statistics_stations %>%
group_by(NODE_NAME, DOY, Parameter) %>%
summarize(Value=mean(Value))
head(Summary_statistics_all)2.3 Plotting
Define order and color scheme for plots
Order for food web plot
ord.pred.fw <- c("Temora", "Centropages", "Pseudocalanus", "Acartia",
"Evadne", "Bosmina", "Synchaeta", "Keratella")
ord.prey.fw <- c("Synechococcales",
"Chlorellales", "Pyramimonadales", # Chlorophytes
"Prymnesiales", # Haptophyta, Hacrobia
"Pyrenomonadales", #Cryptophyta, Hacrobia
"Chromulinales", # Chrysophyceae, Ochrophyta, SAR
"Thalassiosirales", "Chaetocerotales", #Baclillariophyta, Ochrophyta, SAR
"Peridiniales", #Dinoflagellata, SAR
"Nostocales")
color.pred.fw <- c("#8c510a", #Temora
"#bf812d", #Centropages
"#dfc27d", #Pseudocalanus
"#f6e8c3", #Acartia
"#c7eae5", #Evadne
"#80cdc1", #Bosmina
"#35978f", #Synchaeta
"#01665e") #Keratella
color.prey.fw <- c("#3690c0", #Synechococales
"#238443", #Chorellales
"#419753", #Pyramimonadales
"#60AB64", #Prymnesiales
"#7FBF75", #Pyrenomonadales
"#9DD385", #Chromulinales
"#FEA130", #Thalassiosirales
"#FEB23F", #Chaetocerotales
"#ef3b2c", #Peridiniales
"#0570b0" #Nostocales
)Order for line diagrams
ord.pred <- c("Temora", "Centropages", "Pseudocalanus", "Acartia",
"Evadne", "Bosmina", "Synchaeta", "Keratella")
ord.prey <- c("Synechococcales", "Nostocales",
"Thalassiosirales", "Chaetocerotales", #Baclillariophyta, Ochrophyta, SAR
"Peridiniales", #Dinoflagellata, SAR
"Chlorellales", "Pyramimonadales", # Chlorophytes
"Prymnesiales", # Haptophyta, Hacrobia
"Pyrenomonadales", #Cryptophyta, Hacrobia
"Chromulinales" # Chrysophyceae, Ochrophyta, SAR
)
color.pred <- c("#8c510a", #Temora
"#bf812d", #Centropages
"#dfc27d", #Pseudocalanus
"#f6e8c3", #Acartia
"#c7eae5", #Evadne
"#80cdc1", #Bosmina
"#35978f", #Synchaeta
"#01665e") #Keratella
color.prey <- c("#3690c0", #Synechococales
"#0570b0", #Nostocales
"#FEA130", #Thalassiosirales
"#FEB23F", #Chaetocerotales
"#ef3b2c", #Peridiniales
"#238443", #Chorellales
"#419753", #Pyramimonadales
"#60AB64", #Prymnesiales
"#7FBF75", #Pyrenomonadales
"#9DD385", #Chromulinales
"#FEA130", #Thalassiosirales
"#FEB23F", #Chaetocerotales
"#ef3b2c") #Peridiniales
names(color.pred) <- ord.pred
names(color.prey) <- ord.prey
scale_color_pred <- scale_color_manual(breaks = ord.pred,
values = color.pred)
scale_color_prey <- scale_color_manual(breaks = ord.prey,
values = color.prey)
scale_fill_prey <- scale_fill_manual(breaks = ord.prey,
values = color.prey)scale_x_doymonth <- scale_x_continuous('Month',breaks=seq(0,365,by=30.5),
limits=c(0,366),
labels=c("Jan","Feb","Mar","Apr","May","Jun",
"Jul","Aug","Sep","Oct","Nov","Dec"),
expand = c(0,0))Fig. 2 (inner) Annual Fluxes
fw_df <- full_model_df %>%
# Calculate Sum fluxes of all days per year
group_by(Predator, Prey, Station) %>%
summarize(Flux = sum(Flux)) %>%
# Calculate Mean fluxes between station
group_by(Predator, Prey) %>%
summarise(Flux = mean(Flux)) %>%
# Make bipartite matrix format
filter(Flux > 0)
fw_mat <- fw_df %>%
spread(Predator, Flux, fill=0) %>%
# Order the food web with defined orders
relocate(ord.pred.fw) %>% #Columns
arrange(match(Prey, rev(ord.prey.fw))) %>% #Rows
column_to_rownames("Prey") %>%
as.matrix()## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
## # Was:
## data %>% select(ord.pred.fw)
##
## # Now:
## data %>% select(all_of(ord.pred.fw))
##
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.circos.clear()
chordDiagram(fw_mat, grid.col = c(rev(color.prey.fw), color.pred.fw),
annotationTrack = c("grid", "axis"), directional = 1,
direction.type = c("diffHeight", "arrows"), link.arr.type = "big.arrow")
Save plot:
Save summary matrix:
fw_mat %>%
as.data.frame() %>%
rownames_to_column() %>%
write_csv("./Outputs/Tab.S1.csv")Fig. 2 (outer) Fluxes per month
Prepare a food web matrix
fw_df <- full_model_df %>%
mutate(Month = month(as_date(DOY))) %>%
# Calculate Sum fluxes of all days per year
group_by(Predator, Prey, Station, Month) %>%
summarize(Flux = sum(Flux)) %>%
# Calculate Mean fluxes between station
group_by(Predator, Prey, Month) %>%
summarise(Flux = mean(Flux)) %>%
# Make bipartite matrix format
filter(Flux > 0)
head(fw_df)Filter and plot per month
Calculate monthly
fw_df %>%
group_by(Month) %>%
summarize(Flux=sum(Flux)) %>%
mutate(Fluxsqrt=sqrt(Flux))Fig. 3 Selectivity indices
si_mat <- Links %>%
group_by(Predator) %>%
summarize(Preference=Preference/sum(Preference),
Prey = Prey) %>%
mutate(Predator = factor(Predator, ord.pred),
Prey = factor(Prey, ord.prey)) %>%
ggplot() +
geom_bar(aes(x=Predator, y=Preference, fill = Prey), stat = "identity") +
scale_fill_prey +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90))
si_mat
ggsave("./Outputs/fig.3.pdf",plot = si_mat, width = 5, height = 3)Fig. 4 Plot Summary metrics
Fig. 4a Phytoplankton biomass
fig.4a <- Summary_statistics_all %>%
filter(Parameter=="Abundance",
NODE_NAME %in% ord.prey) %>%
ggplot(aes(DOY, Value)) +
geom_line(aes(color=NODE_NAME)) +
scale_x_doymonth +
scale_color_prey +
scale_y_continuous(expression(paste('Biomass(WW) (g / '~m^{2}~')'))) +
labs(color='Prey') +
theme_classic() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
fig.4a
Fig. 4b Zooplankton biomass
fig.4b <- Summary_statistics_all %>%
filter(Parameter=="Abundance",
NODE_NAME %in% ord.pred) %>%
ggplot(aes(DOY, Value)) +
geom_line(aes(color=NODE_NAME)) +
scale_x_doymonth +
scale_color_pred +
scale_y_continuous(expression(paste('Biomass(WW) (g / '~m^{2}~')'))) +
labs(color='Predator') +
theme_classic() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
fig.4b
Fig. 4c Predation losses (Secondary production)
fig.4c <- Summary_statistics_all %>%
filter(Parameter=="Predation_loss",
NODE_NAME %in% ord.prey) %>%
ggplot(aes(DOY, Value)) +
geom_line(aes(color=NODE_NAME)) +
scale_x_doymonth +
scale_color_prey +
scale_y_continuous(expression(paste('Secondary production (kJ / day /'~m^{2}~')'))) +
labs(color='Prey') +
theme_classic() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
fig.4c
Fig. 4d Predation Pressure
fig.4d <- Summary_statistics_all %>%
filter(Parameter=="Predation_pressure",
NODE_NAME %in% ord.prey) %>%
ggplot(aes(DOY, Value)) +
geom_line(aes(color=NODE_NAME)) +
scale_x_doymonth +
scale_color_prey +
scale_y_continuous(expression(paste('Predation pressure (kJ / day / g)'))) +
labs(color='Prey') +
theme_classic() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
fig.4d
Combine and save
fig.4 <- ggpubr::ggarrange(fig.4a, fig.4b, fig.4c, fig.4d)
ggsave("./Outputs/fig.4.pdf",plot = fig.4, width = 8, height = 5)3 Model sensitivit to selectivity
To estimate the importance of zooplankton feeding preferences on annual primary, we resampled the selectivity index data for all predators using the “slice_sample” function in R base package (n=30, resampling with replacement). We also created a null model, assuming no prey preference, and no feeding niche differentiation between the predators, by randomly subsampling the selectivity dataset (n=30). The annual contribution of each primary producer to food web productivity (predation losses) was recalculated with the resampled selectivity indices and was compared with the null model, using linear regression of the log-transformed data.
3.1 Load data, summarize per order
summarized_per_order <- selectivity_data %>%
group_by(NODE_NAME, Order, SAMPLE_date, Sample.x) %>%
summarise(Selectivity_index = sum(Selectivity_index)) %>%
filter(Order %in% ord.prey.fw)## `summarise()` has grouped output by 'NODE_NAME', 'Order', 'SAMPLE_date'. You
## can override using the `.groups` argument.3.2 Subsample and bootstrap selectivity dataset
subsamplePreference <-
function(Subsample) {
names_predator <-
summarized_per_order %>%
spread(Order, Selectivity_index) %>%
pull(NODE_NAME)
names_prey <-
summarized_per_order %>%
spread(NODE_NAME, Selectivity_index) %>%
pull(Order)
subsample <-
summarized_per_order %>%
spread(Order, Selectivity_index) %>%
group_by(NODE_NAME) %>%
slice_sample(prop = 1, replace = TRUE) %>%
gather(key = Prey, value = Preference, Chaetocerotales:Thalassiosirales) %>%
group_by(NODE_NAME, Prey) %>%
summarise(Preference = mean(Preference)) %>%
mutate(Method = "grouped")
random_predator <-
summarized_per_order %>%
spread(Order, Selectivity_index) %>%
ungroup() %>%
slice_sample(prop = 1, replace = TRUE) %>%
mutate(NODE_NAME = names_predator) %>%
gather(key = Prey, value = Preference, Chaetocerotales:Thalassiosirales) %>%
group_by(NODE_NAME, Prey) %>%
summarise(Preference = mean(Preference)) %>%
mutate(Method = "random_predator")
random_prey <-
summarized_per_order %>%
spread(NODE_NAME, Selectivity_index) %>%
ungroup() %>%
slice_sample(prop = 1, replace = TRUE) %>%
mutate(Order = names_prey) %>%
gather(key = NODE_NAME, value = Preference, Acartia:Temora, na.rm = TRUE) %>%
mutate(Prey = Order) %>%
group_by(NODE_NAME, Prey) %>%
summarise(Preference = mean(Preference)) %>%
mutate(Method = "random_prey")
bind_rows(subsample, random_prey) %>%
mutate(Subsample = Subsample)
}3.3 Create extended links dataframe, bootstrapped.
Links_I <- mapply(subsamplePreference, Subsample = 1:30, SIMPLIFY = FALSE) %>%
bind_rows() %>%
mutate(Predator = NODE_NAME) %>%
ungroup() %>%
select(Predator, Prey, Preference, Subsample, Method)
Links_I3.4. Summarise and plot preference bootstrap data
Preference_summary <-
Links_I %>%
group_by(Predator, Prey, Method) %>%
dplyr::summarise(avg_Preference = mean(Preference),
uCI_Preference = Rmisc::CI(Preference)[1],
lCI_Preference = Rmisc::CI(Preference)[3])## `summarise()` has grouped output by 'Predator', 'Prey'. You can override using
## the `.groups` argument.Preference_plot <-
Preference_summary%>%
#filter(Method == "grouped") %>%
mutate(Predator = factor(Predator, ord.pred.fw),
Prey = factor(Prey, ord.prey.fw)) %>%
ggplot(aes(fill=Prey, y=avg_Preference, x=Predator)) +
geom_bar(stat = "identity", position = "dodge")+
geom_errorbar(aes(ymin=lCI_Preference, ymax=uCI_Preference),position = "dodge", size=0.25) +
facet_wrap(Method~.,
labeller = as_labeller(c(`grouped` = "Model",
`random_prey` = "Null-Model"))) +
scale_fill_prey +
theme_bw() +
theme(axis.text.x = element_text(angle = 90)) +
scale_y_continuous(expression(paste('Prey Selectivity Index')))
Preference_plot
3.5 Function: Estimate fluxes of energy
# Wrapper for FoodWebModel function
FoodWebModelI <-
function(nodes=Nodes,
links=Links_I,
biomasses=Biomasses,
temperatures=Temperatures,
doy=152,
station="BY31 LANDSORTSDJ",
iteration=1,
method="grouped") {
links_f <- links %>%
filter(Subsample==iteration,
Method==method)
Model <-
FoodWebModel(nodes = nodes,
links = links_f,
biomasses = biomasses,
temperatures = temperatures,
doy = doy,
station = station) %>%
mutate(Method = method,
Iteration = iteration) %>%
filter(Predator %in% ord.pred.fw,
Prey %in% ord.prey.fw)
return(Model)
}Run the function for all stations, all days, stations, and iterations for the model and the null model.
Method <- c("grouped", "random_prey")
Iteration <- 1:30
Station <- c("BY31 LANDSORTSDJ", "BY15 GOTLANDSDJ", "BY5 BORNHOLMSDJ")
Doy <- 30:340
g <- expand.grid(Method, Iteration, Station, Doy)
# Running this command may be slow, up to 1 h. Concider loading precalculated dataset.
full_model_iteration <-
mapply(FoodWebModelI,
method = g[, 1],
iteration = g[, 2],
station = g[, 3],
doy = g[, 4],
MoreArgs = list(links=Links_I,
biomasses=Biomasses,
temperatures=Temperatures),
SIMPLIFY = FALSE) %>%
bind_rows()
saveRDS(full_model_iteration, "./Outputs/full_model_iteration.rds")full_model_iteration <-
readRDS("./Outputs/full_model_iteration.rds")3.6. Summarize and plot output
summary_full <-
full_model_iteration %>%
mutate(Flux=Flux*86.4) %>% #Modify unit from j/m2/s to Kj/m2/day
group_by(Iteration, Method, Station, Prey) %>%
summarise(Flux = sum(Flux)) %>%
group_by(Iteration, Method, Prey) %>%
summarise(Flux = mean(Flux)) %>%
group_by(Prey, Method) %>%
summarise(avg_Flux = mean(Flux),
uCI_Flux = Rmisc::CI(Flux)[1],
lCI_Flux = Rmisc::CI(Flux)[3]) ## `summarise()` has grouped output by 'Iteration', 'Method', 'Station'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'Iteration', 'Method'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'Prey'. You can override using the
## `.groups` argument.Flux_plot <-
summary_full %>%
mutate(Prey = factor(Prey, ord.prey.fw)) %>%
ggplot(aes(fill=Prey, y=avg_Flux, x=Prey)) +
geom_bar(stat = "identity", position = "dodge")+
geom_errorbar(aes(ymin=lCI_Flux, ymax=uCI_Flux),position = "dodge", size = 0.25) +
facet_wrap(Method~.,
labeller = as_labeller(c(`grouped` = "Model",
`random_prey` = "Null-Model"))) +
scale_fill_prey +
theme_bw() +
theme(axis.text.x = element_text(angle = 90)) +
scale_y_continuous(expression(paste('Energy Flux (kJ / year /'~m^{2}~')')))
Flux_plot
ggpubr::ggarrange(Preference_plot, Flux_plot, ncol = 1, common.legend = TRUE, legend = "right")
ggsave("./Outputs/Fig.S2.pdf")## Saving 7 x 5 in image3.7 Hypothesis tesing - What is the effect of selectivity?
summary_test <-
full_model_iteration %>%
mutate(Flux=Flux*86.4) %>% #Modify unit from j/m2/s to Kj/m2/day
group_by(Iteration, Method, Station, Prey) %>%
summarise(Flux = sum(Flux)) %>%
group_by(Iteration, Method, Prey) %>%
summarise(Flux = mean(Flux)) ## `summarise()` has grouped output by 'Iteration', 'Method', 'Station'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'Iteration', 'Method'. You can override
## using the `.groups` argument.mod <- lm(log(Flux) ~ Prey*Method , data = summary_test)
Anova(mod)plot(mod)
4 Session Info
sessionInfo()## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] car_3.1-0 carData_3.0-5 fluxweb_0.2.0 lubridate_1.8.0
## [5] circlize_0.4.15 svglite_2.1.0 igraph_1.3.5 xts_0.12.2
## [9] zoo_1.8-11 forcats_0.5.2 stringr_1.4.1 dplyr_1.0.10
## [13] purrr_0.3.5 readr_2.1.3 tidyr_1.2.1 tibble_3.1.8
## [17] ggplot2_3.3.6 tidyverse_1.3.2
##
## loaded via a namespace (and not attached):
## [1] fs_1.5.2 bit64_4.0.5 httr_1.4.4
## [4] tools_4.2.1 backports_1.4.1 bslib_0.4.0
## [7] utf8_1.2.2 R6_2.5.1 DBI_1.1.3
## [10] colorspace_2.0-3 withr_2.5.0 gridExtra_2.3
## [13] tidyselect_1.2.0 bit_4.0.4 compiler_4.2.1
## [16] cli_3.4.1 rvest_1.0.3 xml2_1.3.3
## [19] labeling_0.4.2 sass_0.4.2 scales_1.2.1
## [22] systemfonts_1.0.4 digest_0.6.29 rmarkdown_2.17
## [25] pkgconfig_2.0.3 htmltools_0.5.3 dbplyr_2.2.1
## [28] fastmap_1.1.0 highr_0.9 Rmisc_1.5.1
## [31] rlang_1.0.6 GlobalOptions_0.1.2 readxl_1.4.1
## [34] rstudioapi_0.14 shape_1.4.6 jquerylib_0.1.4
## [37] generics_0.1.3 farver_2.1.1 jsonlite_1.8.2
## [40] vroom_1.6.0 googlesheets4_1.0.1 magrittr_2.0.3
## [43] Matrix_1.5-1 Rcpp_1.0.9 munsell_0.5.0
## [46] fansi_1.0.3 abind_1.4-5 lifecycle_1.0.3
## [49] stringi_1.7.8 yaml_2.3.5 plyr_1.8.7
## [52] grid_4.2.1 parallel_4.2.1 crayon_1.5.2
## [55] lattice_0.20-45 haven_2.5.1 cowplot_1.1.1
## [58] hms_1.1.2 knitr_1.40 pillar_1.8.1
## [61] ggpubr_0.4.0 ggsignif_0.6.4 reprex_2.0.2
## [64] glue_1.6.2 evaluate_0.17 modelr_0.1.9
## [67] vctrs_0.4.2 tzdb_0.3.0 cellranger_1.1.0
## [70] gtable_0.3.1 assertthat_0.2.1 cachem_1.0.6
## [73] xfun_0.33 broom_1.0.1 rstatix_0.7.0
## [76] googledrive_2.0.0 gargle_1.2.1 ellipsis_0.3.2